

/*
 * Class:        GammaProcess
 * Description:
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Pierre Tremblay, Jean-Sebastien Parent & Maxime Dion
 * @since        July 2003

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */

package umontreal.iro.lecuyer.stochprocess;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.randvar.*;
import umontreal.iro.lecuyer.probdist.GammaDist;
import umontreal.iro.lecuyer.util.Num;


/**
 * This class represents a <SPAN  CLASS="textit">gamma</SPAN> process
 * 
 * <SPAN CLASS="MATH">{<I>S</I>(<I>t</I>) = <I>G</I>(<I>t</I>;<I>&#956;</I>, <I>&#957;</I>) : <I>t</I>&nbsp;&gt;=&nbsp;0}</SPAN> with mean parameter <SPAN CLASS="MATH"><I>&#956;</I></SPAN> and
 * variance parameter <SPAN CLASS="MATH"><I>&#957;</I></SPAN>.
 * It is a continuous-time process with stationary,
 * independent gamma increments such that for any 
 * <SPAN CLASS="MATH"><I>&#916;t</I> &gt; 0</SPAN>,
 * 
 * <P></P>
 * <DIV ALIGN="CENTER" CLASS="mathdisplay"><A NAME="GammaEqn"></A>
 * <I>S</I>(<I>t</I> + <I>&#916;t</I>) = <I>S</I>(<I>t</I>) + <I>X</I>,
 * </DIV><P></P>
 * where <SPAN CLASS="MATH"><I>X</I></SPAN> is a random variate from the gamma distribution
 * <TT>Gamma</TT>
 * <SPAN CLASS="MATH">(<I>&#956;</I><SUP>2</SUP><I>&#916;t</I>/<I>&#957;</I>, <I>&#956;</I>/<I>&#957;</I>)</SPAN>.
 * 
 * <P>
 * In this class, the gamma process is sampled sequentially using
 * equation.
 * 
 */
public class GammaProcess extends StochasticProcess  {
   // Make sure that process is strictly increasing: when two successive steps
   // are equal, reset the new one to   old*(1 + EPS)
    protected static final double EPS = 1.0e-15;
    protected boolean      usesAnti = false;

    protected RandomStream stream; // used in the generatePath(), so must be kept
    protected GammaGen     Ggen;
    protected double       mu,
                           nu;
    protected double       muOverNu,
                           mu2OverNu;

    // For precomputations for standard GP
    protected double[]     mu2dtOverNu;

   protected void setLarger (double[] path, int left, int mid, int right) {
      /* make sure that path[left] < path[mid] < path[right] by adding a
         small epsilon to the last two */
      double myEps;
      if (path[left] >= 0)
         myEps = EPS;
      else
         myEps = -EPS;
      path[mid] = path[left] * (1.0 + myEps);
      if (path[mid] >= path[right])
         path[right] = path[mid] * (1.0 + myEps);
   }


   protected double setLarger (double[] path, int left, int right) {
      /* make sure that path[left] < s < path[right] by adding a
        small epsilon to the last two; return s. */
      double myEps;
      if (path[left] >= 0)
         myEps = EPS;
      else
         myEps = -EPS;
      double s = path[left] * (1.0 + myEps);
      if (s >= path[right])
         path[right] = s * (1.0 + myEps);
      return s;
   }


   protected double setLarger (double v) {
      /* return a number that is slightly larger than v */
      if (v >= 0)
         return v*(1.0 + EPS) + Num.DBL_EPSILON;
      else
         return v*(1.0 - EPS);
    }


   /**
    * Constructs a new <TT>GammaProcess</TT> with parameters 
    * <SPAN CLASS="MATH"><I>&#956;</I> = <texttt>mu</texttt></SPAN>, 
    * <SPAN CLASS="MATH"><I>&#957;</I> = <texttt>nu</texttt></SPAN> and initial value 
    * <SPAN CLASS="MATH"><I>S</I>(<I>t</I><SUB>0</SUB>) = <texttt>s0</texttt></SPAN>.
    * The gamma variates <SPAN CLASS="MATH"><I>X</I></SPAN> in are generated by inversion
    * using <TT>stream</TT>.
    * 
    */
   public GammaProcess (double s0, double mu, double nu,
                        RandomStream stream) {
        this (s0, mu, nu,  new GammaGen (stream, 1.0));
// Note that the parameters (the 1.0) supplied to
// the GammaGen constructors are meaningless.
// The 'real' parameters are supplied at time of generation
    }


   /**
    * Constructs a new <TT>GammaProcess</TT> with parameters 
    * <SPAN CLASS="MATH"><I>&#956;</I> = <texttt>mu</texttt></SPAN>, 
    * <SPAN CLASS="MATH"><I>&#957;</I> = <texttt>nu</texttt></SPAN> and initial value 
    * <SPAN CLASS="MATH"><I>S</I>(<I>t</I><SUB>0</SUB>) = <texttt>s0</texttt></SPAN>.
    * The gamma variates <SPAN CLASS="MATH"><I>X</I></SPAN> in are supplied by the gamma random
    * variate generator <TT>Ggen</TT>. Note that the parameters of the
    * {@link umontreal.iro.lecuyer.randvar.GammaGen GammaGen} object
    * <TT>Ggen</TT> are not important since the implementation forces the generator
    * to use the correct parameters (as defined above).
    * 
    */
   public GammaProcess (double s0, double mu, double nu, GammaGen Ggen)  {
        this.mu    = mu;
        this.nu    = nu;
        this.x0    = s0;
        this.Ggen     = Ggen;
        stream = Ggen.getStream();
    }

   public double nextObservation()  {
        double s = path[observationIndex];
        double v = s;
        s += Ggen.nextDouble (stream, mu2dtOverNu[observationIndex], muOverNu);
        if (s <= v)
             s = setLarger (v);
        observationIndex++;
        observationCounter = observationIndex;
        path[observationIndex] = s;
        return s;
    } 

   /**
    * Generates and returns the next observation at time 
    * <SPAN CLASS="MATH"><I>t</I><SUB>j+1</SUB> = <texttt>nextTime</texttt></SPAN>,
    * using the previous observation time <SPAN CLASS="MATH"><I>t</I><SUB>j</SUB></SPAN> defined earlier
    * (either by this method or by <TT>setObservationTimes</TT>),
    * as well as the value of the previous observation <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>j</SUB>)</SPAN>.
    * <SPAN  CLASS="textit">Warning</SPAN>: This method will reset the observations time <SPAN CLASS="MATH"><I>t</I><SUB>j+1</SUB></SPAN>
    * for this process to <TT>nextT</TT>. The user must make sure that
    * the <SPAN CLASS="MATH"><I>t</I><SUB>j+1</SUB></SPAN> supplied is 
    * <SPAN CLASS="MATH">&nbsp;&gt;=&nbsp;<I>t</I><SUB>j</SUB></SPAN>.
    * 
    */
   public double nextObservation (double nextT)  {
        // This method is useful for generating variance gamma processes
        double s = path[observationIndex];
        double v = s;
        double previousT = t[observationIndex];
        observationIndex++;
        observationCounter = observationIndex;
        t[observationIndex] = nextT;
        double dt = nextT - previousT;
        s += Ggen.nextDouble (stream, mu2OverNu*dt, muOverNu);
        if (s <= v)
             s = setLarger (v);
        path[observationIndex] = s;
        return s;
    }


   /**
    * Generates, returns and saves the path
    * 
    * <SPAN CLASS="MATH">{<I>X</I>(<I>t</I><SUB>0</SUB>), <I>X</I>(<I>t</I><SUB>1</SUB>),&#8230;, <I>X</I>(<I>t</I><SUB>d</SUB>)}</SPAN>.
    * The gamma variates <SPAN CLASS="MATH"><I>X</I></SPAN> in are generated using the
    * {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream} <TT>stream</TT>
    * or the {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream} included
    * in the {@link umontreal.iro.lecuyer.randvar.GammaGen GammaGen} <TT>Ggen</TT>.
    * 
    */
   public double[] generatePath()   {
        double s = x0;
        double v;
        for (int i = 0; i < d; i++) {
            v = s;
            s += Ggen.nextDouble(stream, mu2dtOverNu[i], muOverNu);
            if (v <= s)
               s = setLarger (v);
            path[i+1] = s;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    } 


   /**
    * Generates, returns and saves the path 
    * <SPAN CLASS="MATH">{<I>X</I>(<I>t</I><SUB>0</SUB>), <I>X</I>(<I>t</I><SUB>1</SUB>),&#8230;, <I>X</I>(<I>t</I><SUB>d</SUB>)}</SPAN>. This method does not use the
    * {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream} <TT>stream</TT> nor the
    * {@link umontreal.iro.lecuyer.randvar.GammaGen GammaGen} <TT>Ggen</TT>. It
    * uses the vector of uniform random numbers <SPAN CLASS="MATH"><I>U</I>(0, 1)</SPAN> provided by the user
    * and generates the path by inversion.  The vector <TT>uniform01</TT> must be of
    * dimension <SPAN CLASS="MATH"><I>d</I></SPAN>.
    * 
    */
   public double[] generatePath (double[] uniform01)  {
        double s = x0;
        double v;
        for (int i = 0; i < d; i++) {
            v = s;
            s += GammaDist.inverseF(mu2dtOverNu[i], muOverNu, 10, uniform01[i]);
            if (v <= s)
                s = setLarger (v);
            path[i+1] = s;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    } 


   /**
    * Sets the parameters  
    * <SPAN CLASS="MATH"><I>S</I>(<I>t</I><SUB>0</SUB>) = <texttt>s0</texttt></SPAN>,
    * 
    * <SPAN CLASS="MATH"><I>&#956;</I> = <texttt>mu</texttt></SPAN> and 
    * <SPAN CLASS="MATH"><I>&#957;</I> = <texttt>nu</texttt></SPAN> of the process.
    * <SPAN  CLASS="textit">Warning</SPAN>: This method will recompute some quantities stored internally,
    * which may be slow if called repeatedly.
    * 
    */
   public void setParams (double s0, double mu, double nu)  {
        this.x0 = s0;
        this.mu = mu;
        this.nu = nu;
        if (observationTimesSet) init(); // Otherwise no need to.
}


   /**
    * Returns the value of the parameter <SPAN CLASS="MATH"><I>&#956;</I></SPAN>.
    * 
    */
   public double getMu()  { return mu; }


   /**
    * Returns the value of the parameter <SPAN CLASS="MATH"><I>&#957;</I></SPAN>.
    * 
    */
   public double getNu()  { return nu; }


   /**
    * Resets the {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream}
    * of the {@link umontreal.iro.lecuyer.randvar.GammaGen GammaGen}
    * to <TT>stream</TT>.
    * 
    */
   public void setStream (RandomStream stream)  {
        this.stream = stream;
        Ggen.setStream (stream);
    }


   /**
    * Returns the {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream} <TT>stream</TT>.
    * 
    */
   public RandomStream getStream()  { return stream; }
 

    protected void init() {
    // init of StochasticProcess only sets path[0] = x0 if observation times are set.
        super.init();
        muOverNu  = mu / nu;
        mu2OverNu = mu * mu / nu;
        mu2dtOverNu = new double[d];
        if (observationTimesSet) {
            for (int i = 0; i < d; i++)
                mu2dtOverNu[i] = mu2OverNu * (t[i+1] - t[i]);
        }
    }

}
